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Abstract 

We propose a fast and efficient algorithm for the detection of metastable states of complex 
dynamical systems from time series. Our approach exploits the characteristic checkerboard texture 
of metastable states exhibited in recurrence plots (RP). On phase space, RPs induce a reflexive 
and symmetric, but not transitive, tolerance relation whose transitive closure yields a phase space 
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partition into metastable states. We construct this partition by a rewriting grammar applied to 
the symbolic dynamics of time indices. 
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The identification of metastable states in complex dynamical systems has become increas- 
ingly important in recent time. Its significance ranges from spin glasses [JJ and molecular 
configurations [2] over geoscientific applications [3] to the neurosciences [IH7]. 

Metastable or quasi- stationary states are regions in the system's phase space with rela- 
tively large dwell that are connected by transients [JJ [8] . Paradigmatic examples are almost 
invariant sets [9], such as the "wings" of the Lorenz attractor centered around its unstable 
foci [10J, or saddles connected by heteroclinic trajectories [71 ITT]. 

For the identification of metastability in time series, their characteristic slow time scales 
have been separated from the fast dynamics of phase space trajectories by several clustering 
algorithms [JJ [2} HJ [91 [TTJ [12]. One method, sometimes called Perron clustering [2], starts 
with some ad hoc partitioning of the system's phase space that leads to an approximate 
Markov chain description [U |2] H] [91 [12] . Applying spectral clustering methods to the 
resulting transition matrix yields the time scales of the process, while their corresponding 
(left-)eigenvectors allow the unification of cells into a partition of metastable states [U IT2"] . 
Another approach by Hutt et al. [UJ utilizes the slowing-down of the system's trajectory in 
the vicinity of saddles by means of phase space clustering. 

However, those methods might be computationally cumbersome in some applications. In 
particular, spectral clustering and partitioning could be problematic for higher-dimensional 
spaces and sparse data, because either the full phase space or a lower-dimensional projection 
has to be initially partitioned into equally populated cells from which counting measures are 
to be derived for the estimation of Markov transition matrices [I] . 

In this Letter, we propose a parsimonious, fast and feasible algorithm for the identification 
of metastable states in dynamical systems from measured or simulated time series. Our 
approach is inspired by Poincare's famous recurrence theorem [13]: When B £ (xq) is a "ball" 
of radius e > centered at initial condition xq G X, where X C M d denotes the system's 
phase space of dimension d, then there is an infinite sequence of iterates x tk = & k (xo) 
returning to B e (xq), such that Xt k G B e (xo) for all sub- indices fceN. Here, $* : X — > X for 
a given time t describes the system's flow in phase space. These recurrences can be visualized 
by means of Eckmann et al.'s [14J recurrence plot (RP) method where the element 



of the recurrence matrix R = (Rij) is one if Xj G B £ (xi) and zero otherwise [TH [15], as 
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mediated by the Heaviside step function 0. Eckmann et al. [H] have already pointed out 
that RPs display metastability by a characteristic "checkerboard texture". We illustrate 
this in Fig. [T] with two paradigmatic RPs for the Lorenz attractor [TU] and for a stable 
heteroclinic contour of three competing Lotka-Volterra populations [7] in Fig. [Tj^a) and (b), 
respectively [2T] . 
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Figure 1: Recurrence plots [Eq. Q] for Euclidian norm of (a) Lorenz attractor and (b) stable 
heteroclinic contour of a Lotka-Volterra system. Parameters: (a) e = 5, (b) e = 0.1. Black pixels 
denote Rij = 1, white ones Rij = 0. 

Both RPs exhibit the typical texture of diagonal line patters that are characteristic for 
oscillatory dynamics. For the Lorenz system shown in Fig. [T](a) these oscillators correspond 
to the metastable wings. Going along the line of identity (LOI) reveals transient transitions 
between about four metastable states, while the horizontal and vertical recurrences of the 
same patterns indicate that there are indeed only two metastable states being involved, 
namely the wings, repeatedly explored by the system's trajectory. On the other hand, 
the texture of Fig. [TJb) represents a periodic winnerless competition [16] where the 
characteristic slowing-down in the vicinity of a saddle [TT] is reflected by the broadening of 
horizontal, respective, vertical lines, indicating an increased number of recurrences within 
an e-ball around the saddle [T5] . 

For uniform e, the recurrence matrices obtained from Eq. are reflexive, Ra = 1 (the 
LOI), and symmetric, = Rji, but in general not transitive, i.e. R^ = 1 and Rjk = 1 do 
not necessarily imply R ik = 1. In order to cope with this disadvantage, Faure and Lesne 
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[T7] suggested to compute the recurrence matrix from words in a symbolic dynamics [18] 
through 

Rij = &WiWj 1 (2) 

where Wi, Wj are words of length m at times i and j in a symbolic sequence s = a\a^ ■ ■ ■ a n . 
Here, 5 a & = 1 if a = b and zero otherwise, denotes the Kronecker matrix. Symbolic RPs 
given by Eq. ^ are also transitive, because symbolic dynamics results from a partition of 
the system's phase space into equivalence classes. 

However, for the more general case of e-RPs [Eq. Q], we can still define a tolerance 
relation xRy [TH] on phase space X as follows: we say x G X i?-tolerates y G X, xRy, if 
there exists an initial condition xq G X and a time t such that (1) x G B £ (xq), (2) y G B £ (xt) 
and (3) x t G B e (xo), i.e. i?ot = 1- Here, x t = &(xo) is again the iterate of xq. The tolerance 
relation R is reflexive: xRx, since i?oo = 1; it is symmetric: xi??/ implies yRx as i?o* = Rto] 
but it is not transitive: xRy and do not imply xRz 

Nevertheless, we can construct the transitive closure, R + , of R which is possible for 
every relation. This new relation would be an equivalence relation whose equivalence classes 
partition the phase space into unions of recurrent and thus intersecting e-balls which cover 
the metastable states. We achieve this construction by regarding the e-RP R [Eq. ([T])] as a 
grammatical rewriting system over the time indices of a given trajectory xt [20]. Thus, we 
first map the trajectory x t to the sequence of successive time indices, regarded as symbols: 
Xt — > s t = t. Then we define a formal grammar of rewriting rules: if i > j and Rij = 1 create 
a rule % — > j. To enforce transitivity for i > j > k, R^ = 1 and Rn~ = 1, we first eliminate 
the redundancy by rewriting only % — > k and then create an additional rule j — > k. Finally, 
we apply this grammar to the initial sequence of time indices St = t in order to replace large 
indices by smaller ones, thus exploiting the recurrence structure of the data. The result is a 
transformed symbolic sequence s' t , whose symbolic RP R + [Eq. [T7] yields the transitive 
closure R + of the tolerance relation R on phase space. 

As R + is an equivalence relation, its equivalence classes partition the phase space into 
metastable states. This partition contains the union of intersecting £-balls along the sampled 
trajectory together with its complements. 

Let us illustrate the procedure by means of a simple example. Assume, we have a series 
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of only five data points (xi, X2, ■ ■ ■ , X5) that gave rise to the recurrence matrix 
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The algorithm starts in the 5th row, detecting a recurrence R52 = 1- Since 5 > 2, we create a 
rewriting rule 5 — > 2. Because the next recurrence in row 5 is trivial, the algorithm continues 
with row 4, where R41 = i? 42 = 1. Now, two rules 4 — > 1 and 4 — > 2 could be generated. 
However, the latter is redundant. Therefore the algorithm only records the rule 4 — ^ 1 . 
Moreover, transitivity is taken into account by an additional rule 2 — > 1. Next, row 3 does 
not contribute to the algorithm and rows 2 and 1 can be neglected due to the symmetry. 
Recursively applying this grammar to the symbolically encoded time series s = <cl2345» 
yields s' = «11311», i.e. a bistable system with states «1» and «3». 

In order to validate our construction, we employ the method to the examples from Fig. [TJ 
Yet, for the ball size e chosen in Fig. [TJ every phase space point turns out to be equivalent to 
every other point, leading to a trivial partition. We therefore diminish e in such a way that 
the resulting e-RPs become too sparse for visualization. Furthermore, we find a tradeoff 
between e and the sampling rate. When the latter is too small, metastable states could be 
reflected by spurious oscillations in the symbolic dynamics s' t . Thus, we also increased the 
sampling rate for the Lorenz attractor and the Lotka-Volterra system. 

Figure [2] shows the three time series X2(t), Xs(t) from the numeric solution of the 

Lorenz attractor (a) and of the heteroclinic contour in the Lotka-Volterra populations (b), 
together with the suitably rescaled symbolic sequence s' obtained from our recurrence-based 
partitioning method. 
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Figure 2: (Color online) Time series of (a) Lorenz attractor and (b) stable heteroclinic contour 
of a Lotka-Volterra system. Blue: x\(t), green: X2(t), red: xs(t). The black curve displays the 
(rescaled) symbolic sequence s'(t) from the recurrence-based phase space partition. Parameters: 
(a) e = 1, (b) e = 0.013. 

In Fig. the switching from one metastable state, i.e., wing, into another one is 

simultaneously reflected by the sign change of coordinates X\, x 2 and by the jumps between 
the long-lasting plateaus of the symbolic dynamics s'. Short transients between metastable 
states are represented by monotonicaliy increasing pieces of s'. This is even more salient 
in Fig. ^b), where most of the heteroclinic dynamics is transient. On the other hand, the 
characteristic slowing down in the vicinity of the saddles leads to increased recurrences and 
hence to the plateaus in s', again. 

Figure [3] displays the symbolic recurrence plots Eq. ([2J of s' for word length m = 1, -R + , 
as the transitive closure R + of the tolerance relation R. 
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Figure 3: Symbolic recurrence plots [Eq. ([2])] of (a) Lorenz attractor and (b) stable heteroclinic 
contour of a Lotka-Volterra system. Black pixels denote Rfj = 1, white ones Rfj = 0. 

After the recurrence-based symbolic encoding, the checkerboard texture indicating the 
metastable states is highly enhanced for the Lorenz attractor Fig. |3](a) in comparison to 
Fig. [TJa) and for the Lotka-Volterra system Fig. ^b) compared with Fig. [ijb). 

Finally, we combine the phase space information contained in the state vector x G X 
with that contained in the symbolic dynamics s' at one instance of time t for identifying 
metastable states in phase space. 

Figure [4] depicts the trajectories around the Lorenz attractor (a) and around the hetero- 
clinic contour (b) where each sampled state x t at time t is colored according to the symbol 
s' t using a rotating color palette. 
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Figure 4: (Color online) Phase space partitions into metastable states obtained from recurrence- 
based symbolic dynamics, (a) Lorenz attractor and (b) stable heteroclinic contour of a Lotka- 
Volterra system. Colored balls indicate the membership to a partition cell. 

Now, transients are represented by continuously changing colors while metastable states 
are homogeneously colored regions in phase space. For the Lorenz attractor in Fig. |4^a) 
these are the wings essentially colored in turquoise and dark blue. For the Lotka-Volterra 
dynamics in Fig. [4](b) the metastable states are the vertices of the "triangle" colored in 
blue, red and purple. 

In this Letter we proposed a parsimonious, fast and feasible algorithm for the detec- 
tion of metastable states in complex dynamical systems. In contrast to techniques based 
on Markov chains, which require an ad hoc partitioning of the system's phase space into 
equally populated cells, the estimation of transition probabilities by counting measures and 
subsequent spectral clustering methods, our approach simply exploits the recurrence struc- 
ture of the system's dynamics by taking the transitive closure of the recurrence plot, thereby 
partitioning the phase space into unions of intersecting e-balls along the actual trajectory. 

The proposed method could have a number of interesting applications in many different 
fields, such as molecular dynamics [2], geo- [3], and neurosciences jlHTj for the identification 
of metastable behavior. Moreover, it could also be useful for the analysis of complex networks 
for solving graph partition and related problems by taking the transitive closure of the 
graph's adjacency matrix. 
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